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Abstract 

Using the coupled cluster method we study the phase diagram of the spin- 1/2 Heisenberg antifer- 
romagnet on a honeycomb lattice with nearest-neighbor exchange coupling J\ > and frustrating 
next-nearest-neighbor coupling J2 = xJ\ > 0. In the range < x < 1 we find four phases exhibiting 
respectively Neel, plaquette, staggered dimer, and Neel-II orderings, with quantum critical points 
at x Cl ph 0.207(3), x C2 sa 0.385(10), and x c . A pa 0.65(5). The transitions at x Cl and x C3 appear to be 
continuous (and hence deconfined) ones, while that at x C2 appears to be a direct first-order one. 

PACS numbers: 75.10.Jm, 75.10.Kt, 75.30.Kz, 75.50.Ee 
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Frustrated quantum spin models on the two-dimensional (2D) honeycomb lattice have 
become the objects of intense study. Impetus came from the reported presence of a quantum 
spin-liquid (QSL) phase in both the exactly soluble (albeit somewhat artificial) Kitaev model 
of spin-1/2 particles on a honeycomb lattice [lj], and the half-filled Fermi-Hubbard (FH) 
model on a honeycomb lattice [2j. Thus, Meng et al. reported in a quantum Monte Carlo 
(QMC) calculation, free of the usual fermion sign problems, the presence in the honeycomb 
FH model of a QSL phase, at moderate values of the on-site Coulomb repulsion strength 
(U), situated between the nonmagnetic metallic insulator (or semi- metal) phase at low U 
and the antiferromagnetic (AFM) Mott insulator phase for large U. Since the U — > oo limit 
corresponds to the pure Heisenberg antiferromagnet (HAFM), i.e., with nearest-neighbor 
(NN) interactions (of strength J\ > 0) only, the Mott insulator phase of the Hubbard 
model corresponds to the Neel-ordered phase of the HAFM spin-lattice model. Higher-order 
terms in the t/U expansion of the FH model (where t is the strength parameter of the 
NN hopping term) lead to frustrating exchange couplings in the corresponding spin-lattice 
model in which the HAFM with NN exchange couplings is the leading term in the large- £/ 
expansion. The simplest such frustrated model is the J\—Ji model studied here, where the 
next-nearest-neighbor (NNN) spin pairs have an additional exchange coupling of strength 
J 2 > 0. 

A later study of the FH model, using a Schwinger boson mean field theory (SB-MFT) 
approach Q], provided some corroboratingevidence for a Z 2 QSL state; and a Schwinger 
fermion representation of the same model [4j gave some evidence for both a Z 2 QSL phase 



and a chiral anti 
by Sorella et al. 



'erromagnetic phase. However, later numerically exact QMC calculations 



5|, with much larger clusters than those used by Meng et al. j^j], have cast 
considerable doubt on their original finding of an intermediate QSL phase. Nevertheless, 
the study of analogous frustrated spin models on the honeycomb lattice remains of great 
interest in its own right. Thus, at intermediate values of U, around those for which Meng 

n 

et al. [2[ initially reported the QSL, the properties of the FH model are mirrored, at least 
loosely, by a spin-1/2 J\—Ji model on the honeycomb lattice. For this and other reasons, 
this model and its generalizations [specifically to include also next-next-nearest-neighbor 
(NNNN) bonds with strength J 3 ], have been much studied recently. Its Hamiltonian 
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FIG. 1: (Color online) The J1-J2 model on the honeycomb lattice (with J\ = 1), showing (a) the 
Neel and (b) Neel-II states. The arrows represent spins located on lattice sites •. 



(1) 



is given by 

H = Ji } j Si ■ Sj + J 2 ^ s.j • s fc , 

where index i runs over all honeycomb lattice sites, and indices j and k run over all NN 
and NNN sites to i, respectively, counting each bond once only. Each lattice site i carries 
a particle with spin s = \ and a spin operator Sj = (sf,sf,sf). The lattice and exchange 
bonds are illustrated in Fig. [TJ We are interested in the case where both NN and NNN 
bonds are AFM in nature, and henceforth, we put J\ — 1 to set the energy scale and define 
the frustration parameter x = J2/ J\- 

The classical (s — > 00) ground state of the model is Neel-ordered for < x < g, whereas 
for all values x > | the spins are spirally ordered. In this latter regime, the classical 
model has a one-parameter family of degenerate incommensurate ground states where the 
spiral wave vector can orient in any direction. At leading order, i.e., 0(1/ s), spin-wave 
fluctuations lift this accidental degeneracy in favor of particular wave vectors p. For the 
extreme quantum case, s = 1/2, considered here, we expect quantum fluctuations to be 
strong enough to destroy the spiral order over a wide range of values of x. In a recent paper 



12] that used the coupled cluster method (CCM), we have verified that expectation for all 



values in the range < x < 1 considered here. 



We showed too 



12| that quantum fluctuations preserve the Neel order to higher values 



of x than in the classical model. Thus, we found that the ground-state (GS) phase of the 
s = 1/2 model is Neel-ordered for x < x Cl ~ 0.207(3). At x — x Cl there appears to be 
a continuous deconfmed phase transition to a GS paramagnetic phase exhibiting plaquette 
valence-bond crystalline (PVBC) order. Furthermore, we found the PVBC state to be the 
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stable GS phase in the regime x Cl < x < x C2 , where x C2 ~ 0.385(10). Our aim now is to 
investigate further the transition at x — x C2 and the nature of the GS phase(s) for x > x C2 . 

The CCM, that we will employ here, has been very successfully applied to many models 
in quantum magnetism, including models on the honeycomb lattice 131 1 of interest 

here. It provides a well-structured means of studying various candidate GS phases and their 
regimes of stability, for each of which the description is systematically improvable in terms of 
well-defined truncation hierarchies for the quantum multi-spin correlations. We now briefly 
describe the method and refer the reader to Refs. 
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for further details. 



The starting point for any CCM calculation is the selection of a suitable normalized 
model (or reference) state |$). For spin systems it is often convenient to take a classical 
(uncorrelated) GS wave function for |$). For the present case we choose the Neel state 
shown in Fig. []Ta) for small values of the frustration parameter x. For larger values of x 
we could choose one of the classical spiral GS phases to provide a CCM model state, but 
as we have argued above these are likely to be very fragile against quantum fluctuations. 
Instead, for larger values of x, we choose here the so-called Neel-II phase shown in Fig. Hlb) 
(which has also been denoted as the anti-Neel phase earlier 12[), that occurs in the classical 
(s — > oo) model only at the isolated and highly degenerate critical point x = \. Whereas 
the Neel state has all 3 NN spins to a given spin antiparallel to it, the Neel-II state also 
comprises AFM sawtooth chains along one of the three equivalent honeycomb directions, 
but with NN spins on adjacent chains now parallel to one another. The Neel-II state thus 
breaks the lattice rotational symmetry. 

It is convenient to perform a mathematical rotation of the local axes of the spins such that 
all spins in the reference state align along the negative z-axis. The Schrodinger ground-state 
ket and bra CCM equations are H\^l) = E\ty and (&\H = E$\ respectively. The CCM 
employs the exponential parametrizations, = e 5 |$) and = ($|iSe~ 5 . The correlation 
operator S is expressed as S = JZi^o^i^f an d its counterpart is S — l + X^o^^F where, 
by definition, Cf\$) = = ($1(7/, V/ ^ 0. Thus we have the normalization condition 
(^1^) = ($1$) = 1- The multispin creation operators Cj = {Cj) with Cq = 1, are 
written as Cf = s^st . . ■ s^, in terms of the single-site spin-raising operators s£ = sl + is y k . 
The GS energy is E = ($|e _5 ife 5 |$); and the local average onsite magnetization M in 
the rotated spin coordinates is M = — ^(^| Ylf=i s j\^)- The ket- and bra-state correlation 
coefficients (5/, Si) are calculated by requiring the expectation value H = (&\H\ty) to be 
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a minimum with respect to all parameters (Si, Si), and hence (3>|Cj e s He s \Q) = and 
($\S(e- s He s - E )C+\<S>) = ; V/ + 0. 

The CCM formalism is exact if all spin configurations are included in the S and S op- 
erators. In practice, however, truncations are needed. We em ploy here the well-studied 
localized (lattice-animal-based subsystem) LSUBm scheme 18M23||. in which all possible 
multi-spin-flip correlations over different locales on the lattice defined by m or fewer con- 
tiguous lattice sites are retained. Such clusters are defined to be contiguous in this sense 
if every site in the cluster is adjacent (as a nearest neighbor) to at least one other site in 
the cluster. The numbers Nf of such fundamental configurations that are distinct under the 
(space and point-group) symmetries of the lattice and the model state increase rapidly with 
the LSUBm truncation index m. Thus the highest LSUBm level that we can reach here, 



24|, is LSUB12, 



even with massive parallelization and the use of supercomputing resources 
for which N f = 293309 for the Neel-II state. 

Since, in any truncation, CCM parametrizations automatically satisfy the Goldstone 
linked cluster theorem, we may work from the outset in the thermodynamic limit, N — > oo. 
Nevertheless, the raw LSUBm data still need to be extrapolated to the exact m — > oo limit. 
Although there are no known exact extrapolation rules we have much experience in doing 
so. Thus, for the GS energy per spin, E/N, we use (see, e.g., Refs. 
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2&) 



E(m)/N = ao + aim 2 + 02m 



(2) 



while for the magnetic order parameter, M, defined above, we use either the scheme 

M(m) = b + bxm' 1 + b 2 m~ 2 , 



(3) 



for systems showing no or only slight frustration (see, e.g., Refs. 19j, |2l|), or the scheme 



M(m) = c + cim~ 1/2 + b 2 m~ 3/2 



(4) 



for more strongly frustrated systems or ones showing a GS order- disorder transition (see, 



e.g., Refs. 22|, |23j). Since the hexagon is an important structural element of the honeycomb 



lattice we perform the extrapolations here only for LSUBm data with m > 6. 

In Fig. [2] we show our results for the magnetic order parameter, M, using both the 
Neel and Neel-II states as CCM model states. On the Neel side we show extrapolations 
using both Eqs. Q and Equation Q, which is appropriate when J 2 = 0, yields 
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FIG. 2: (Color online) CCM LSUBm results for the GS order parameter, M, of the spin-1/2 J1-J2 
honeycomb model (J\ = 1), using the Neel (left curves) and Neel-II (right curves) states as model 
states, with m = {6,8, 10, 12}. The extrapolated curves LSUBoo(l) and LSUBoo(3) use this data 
set with Eqs. (jU) and (|3|) respectively, while the LSUBoo(2) curve uses Eq. ([3]) with the restricted 
set m = {8, 10,12}. 

the value M ~ 0.271(2) for the unfrustrated HAFM on the hexagonal lattice (i.e., with 
NN interactions only), in excellent agreement with the best available QMC estimate j^J, 
M = 0.2677(6). Our own error estimates are based on sensitivity checks using different 
LSUBm data sets. Similarly we see from Fig. [2] that all extrapolations give essentially the 
same estimate x Cl ~ 0.207(3) for the point where Neel order vanishes (M — > 0). We showed 



previously [12| that the phase transition at X — X IS Ci continuous deconfmed one between 
states with Neel and PVBC order. 

Figure [2] also shows corresponding results for M for a possible phase with Neel-II ordering. 
Inspection shows clearly that inclusion of the LSUBm result with m = 6 for this state is not 
compatible with the expected extrapolation scheme of Eq. (j3J), whereas those with m > 6 
fit very well. Why the m = 6 result should be anomalous in this case is unclear. For these 
reasons we show in Fig. [2] only extrapolated results using the Neel-II model state based on 
m = {8, 10, 12}. The results clearly show that Neel-II ordering is present, albeit with a 
rather small value of the order parameter, M < 0.1, for x > x Cz where x Cz ~ 0.65(5), but 
where the error estimate is now very rough and rather uncertain. 

In our previous work [3| we showed that the Neel-II state becomes susceptible to PVBC 
ordering for x < x C2 ~ 0.385(10), but we now observe that the Neel-II state is itself only 
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FIG. 3: (Color online) Left: CCM LSUBm results for the inverse staggered dimer susceptibility, 
of the spin-1/2 J\-Ji honeycomb model (Ji = 1), using the Neel-II state as model state, with 
m = {6, 8, 10, 12}. The extrapolated curves LSUBoo(l) and LSUBoo(2) are derived from fitting the 
perturbed energies (see text) as e(5) = eo(5) + e\{8)m~ v , and use the data sets m = {6,8, 10, 12} 
and m = {8, 10, 12} respectively. Right: The field F — > 5 Od for the staggered dimer susceptibility, 
Xd- Thick (red) and thin (black) lines correspond respectively to strengthened and unaltered NN 
exchange couplings, where Od = j) a ij s i' s j, an d the sum runs over all NN bonds, with a%j = +1 
and for thick (red) lines and thin (black) lines respectively. 
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stable as a magnetically ordered state for x > x C3 . We are thus led to enquire about the 
possible GS phase(s) of the system in the range x C2 < x < x C:j . In view of the persistence 
of our CCM LSUBm solutions based on the Neel-II model state, with finite values of m, 
well into the region x < x C3 before they terminate (as is clearly seen from Fig. |2]), we expect 
that the actual GS phase in this intermediate regime might share similarities with the Neel- 
II state. For example, just as the Neel-II state breaks the lattice rotational symmetry, so 
does another valence-bond solid state, namely the staggered- dimer valence-bond crystalline 
(SDVBC) (or lattice nematic) state. This is formed from the Neel-II state by replacing all 
of the parallel NN spin pairs by spin-zero dimers (and see Fig. [3]). 

In order to investigate the possibility of an SDVBC phase we first consider the response 
of the system to a field operator F (and see, e.g., Ref. |23[). Thus, a field term F = 5 Od is 
added to the Hamiltonian of Eq. <^f, where Od is an operator corresponding to the possible 
SDVBC order, illustrated in Fig. [3]and denned in its caption. The energy per site, E(S)/N = 
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e(S), is then calculated in the CCM for the perturbed Hamiltonian H + F, using the Neel- 
II model state. We define the corresponding susceptibility as Xd = — {d 2 e(S)) / (d5 2 )\ s=0 . 
Clearly the GS phase becomes unstable against SDVBC order when Xd~ l becomes zero. We 
now use the LSUBm extrapolation scheme e(5) = eo(5) + ei(5)m~' / , with the exponent v 
also a fitting parameter, rather than our standard energy extrapolation scheme of Eq. 
in order to calculate the extrapolated values of x~d 1 shown in Fig. [3j For the same data set 
m = {8, 10, 12} used to calculate M for the Neel-II state above, the fitted value of v is close 
to 2 over most of the range of the J2 values shown, except near the termination point of 
this phase, where it falls sharply. By contrast, for the set m = {6, 8, 10, 12} also shown in 
Fig. [31 v is closer to 1 over most of the range. This again reinforces the anomalous nature 
of the LSUB6 results. 

What we see from Fig. [3] is that the extrapolated value of Xd 1 * s c l° se to zero over a 
range of values of x that extends from x C2 below to an upper value of about 0.6, which 
is completely compatible with the value x C3 obtained from the order parameter M of the 
Neel-II state. Thus, by combining our results, we conclude that in the region x C2 < x < x cz 
the GS phase has SDVBC order, while for x > x C3 the GS phase has Neel-II order, although 
this latter ordering is weak and quite fragile against the still strongly competing SDVBC 
order. The shape of the CCM curves for x^" 1 in Fig. [3] are indicative of a continuous (and 
hence deconfined) quantum critical point at x C:j , whereas the corresponding curves for x" 1 , 
the inverse plaquette susceptibility, found in our earlier work [l2| were much more indicative 
of a direct first-order transition at x C2 . We see no signals at all of the spiral ordering that is 
present classically for x > g for any value of x in the range < x < 1 examined. 

In conclusion, over the range < x < 1 we find that the spin-1/2 J\—Ji HAFM on the 
honeycomb lattice has four phases with, respectively, Neel, PVBC, SDVBC, and Neel-II 
ordering. Our first calculated critical point, x Cl ~ 0.207(3), at which Neel order melts, 
agrees well with other recent results, including x ci ~ 0.195(25) from a large-scale exact 



diagonalization (ED) study |8j, x ci ~ 0.26 [14j and x Cl ~ 0.22 [15| from two separate 



density-matrix renormalization group (DMRG) studies, and x Cl ~ 0.2075 from a SB-MFT 



approach 



16|. Both DMRG studies 14|, ll5| and the ED study [8| concur with us that the 



transition at x r _, is probably a continuous deconfined one to a PVBC state, whereas the SB- 



MFT study 



16] indicates a transition to a QSL state. Our second calculated critical point, 



x C2 ~ 0.385(10), at which the PVBC order melts, is similarly in good agreement with the 



S 



result x C2 ~ 



0.375(25) from the ED study |8|, and the results x 



<:■> 



0.36 



and x, 



C-2 



0.35 



15[ from the two DMRG studies. We find that the transition at x r „ is probably a direct 
first-order one to a state with SDVBC order. Both DMRG studies 14. 1 15 1 concur that the 
transition at x C2 is to a state with SDVBC order, although Ganesh et al. 15| find evidence 
for the surprising scenario that the transition at x C2 is also of the continuous deconfined type, 
as at x Cl . The ED study 8|, by contrast, finds a first-order transition at x C2 to a state that 
cannot be distinguished between having either SDVBC or Neel-II order. Finally, we find 
evidence for a third critical point at x C:j ~ 0.65(5) at which a continuous (and hence again 
deconfined) transition occurs to a state with weak Neel-II magnetic order. It is interesting 
to note that such a transition is also compatible with the DMRG result of Ganesh el al. 



which could not rule out a melting of the SDVBC order for values x > 0.7. 
We thank the University of Minnesota Supercomputing Institute for the grant of super- 
computing facilities for this research. 
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